function out = getEquityPrem(model,modelRealYields,setupFilter)

outSim  = simProjection(model,setupFilter.shocks);
posMuz  = find(strcmp(model.labelx,'$\mu _{z,t}$')); 
posDiv  = find(strcmp(model.labely,'$div_t$')); 
posEq   = find(strcmp(model.labely,'$p^m_t$')); 
logmuz_cu = outSim.x_cu(posMuz,:)+log(model.params.MUZss);
div_cu  = outSim.y_cu(posDiv,:);
eq_cu   = outSim.y_cu(posEq,:);
% Equity return 
retEq_cu= log([1 (exp(div_cu(2:end))+exp(eq_cu(2:end)))./exp(eq_cu(1:end-1))])+logmuz_cu;
% Real yield
if model.pruningOn == 1
    realYields_cu = getYields(modelRealYields,outSim.xAll_cu);
    yields_cu     = getYields(model,outSim.xAll_cu);
elseif model.pruningOn == 0
    realYields_cu = getYields(modelRealYields,outSim.x_cu);
    yields_cu     = getYields(model,outSim.x_cu);
end

% Equity premium
levFactor   = 2;
exReturn_cu = levFactor*(retEq_cu-realYields_cu(1,:));
out.meanEqPrem  = mean(exReturn_cu*4);
out.stdEqPrem   = std(exReturn_cu)*sqrt(4);
out.meanRetEq   = mean(levFactor*retEq_cu*4);
out.stdRetEq    = std(levFactor*retEq_cu)*sqrt(4);

% Computing annual price-dividend ratio
numSim          = size(exReturn_cu,2);
logz_cu         = cumsum(logmuz_cu-log(model.params.MUZss));
% Cannot include the deterministic trend because we then get overflow
% The leverage factor is 2 i.e. div^2
div_quarterly   = exp(div_cu).*exp(logz_cu);
div_annual      = zeros(1,numSim);
exReturn_annual = zeros(1,numSim);
for s=1:numSim
    if s == 1
        div_annual(1,s)      = div_quarterly(1,s)*4;
        exReturn_annual(1,s) = exReturn_cu(1,s)*4;
    elseif s == 2
        div_annual(1,s)      = (div_quarterly(1,s)+div_quarterly(1,s-1))*2;
        exReturn_annual(1,s)   = (exReturn_cu(1,s)+exReturn_cu(1,s-1))*2;
    elseif s == 3
        div_annual(1,s) = div_quarterly(1,s)+div_quarterly(1,s-1)...
            +div_quarterly(1,s-2)*4/3;
        exReturn_annual(1,s) = (exReturn_cu(1,s)+exReturn_cu(1,s-1)+exReturn_cu(1,s-2))*4/3;
    elseif s>3
        div_annual(1,s) = div_quarterly(1,s)+div_quarterly(1,s-1)...
            +div_quarterly(1,s-2)+div_quarterly(1,s-3);
        exReturn_annual(1,s) = (exReturn_cu(1,s)+exReturn_cu(1,s-1)+exReturn_cu(1,s-2)+exReturn_cu(1,s-3));
    end
end
pd              = exp(eq_cu).*exp(logz_cu)./div_annual;    %the deterministic trend cancel approximately out
out.meanLogPd   = mean(log(pd));
out.stdLogPd    = std(log(pd));

% Predicting excess returns by the slope of the yield curve (and pd)
maxHorizon         = 12;
regExR_on_pd       = zeros(maxHorizon,2);
regExR_on_slope    = zeros(maxHorizon,2);
regExR_on_slope_pd = zeros(maxHorizon,3);
%exR_z              = (exReturn_cu'-mean(exReturn_cu))./std(exReturn_cu);
exR_z              = (exReturn_annual'-mean(exReturn_annual))./std(exReturn_annual);
log_pd_z           = (log(pd)'-mean(log(pd)))/std(log(pd));
slope              = (yields_cu(end,:)-yields_cu(4,:))*4;
slope_z            = (slope'-mean(slope))/std(slope);
for j=1:maxHorizon
    results = predictRegress(exR_z,log_pd_z,j,j*2);
    regExR_on_pd(j,:)   = [results.beta(2,1) results.rsqr];
    results = predictRegress(exR_z,slope_z,j,j*2);
    regExR_on_slope(j,:)   = [results.beta(2,1) results.rsqr];
    results = predictRegress(exR_z,[slope_z log_pd_z],j,j*2);
    regExR_on_slope_pd(j,:)   = [results.beta(2,1) results.beta(3,1) results.rsqr];
end
%select = [4,8,12]';
%regExR_on_slope(select,:)
%regExR_on_slope_pd(select,:)
out.regExR_on_pd = regExR_on_pd;
out.regExR_on_slope = regExR_on_slope;
out.regExR_on_slope_pd = regExR_on_slope_pd;


end